Script in series: “Associations between early life mental health indicators and sleep in adulthood”
This script performs multiple imputation and regression analyses on the dataset associated with this project.
We start by loading the necessary dependencies.
library(data.table)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.1 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between() masks data.table::between()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks data.table::first()
## ✖ lubridate::hour() masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::last() masks data.table::last()
## ✖ lubridate::mday() masks data.table::mday()
## ✖ lubridate::minute() masks data.table::minute()
## ✖ lubridate::month() masks data.table::month()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ lubridate::second() masks data.table::second()
## ✖ purrr::transpose() masks data.table::transpose()
## ✖ lubridate::wday() masks data.table::wday()
## ✖ lubridate::week() masks data.table::week()
## ✖ lubridate::yday() masks data.table::yday()
## ✖ lubridate::year() masks data.table::year()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(magrittr)
##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extract
library(parallel)
library(mice)
##
## Attaching package: 'mice'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## cbind, rbind
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(sandwich)
library(nnet)
We now load the cleaned BCS70 data and prepare them for multiple imputation.
data_merged <- fread(".\\Cleaned\\data_clean.csv") %>% data.frame()
# Replacing blank values in character vectors with NA
is.na(data_merged[,]) <- data_merged[,] == ""
# Linearly rescaling all numeric variables to between 0 and 1 to aid regression model convergence
data_merged %<>% mutate_if(is.numeric, function(x) (x - min(x, na.rm = TRUE))/(max(x, na.rm = TRUE) - min(x, na.rm = TRUE)))
# Listing dependent variables for Poisson and multinomial regression, respectively
pois_list <- c("dv_sr_sleep_abn", "dv_sd_sleep_abn", "dv_pal_sleep_abn", "dv_vdb_sleep_abn", "dv_winkler_sleep_abn")
multinom_list <- c("dv_sr_sleep_abn_cat", "dv_sd_sleep_abn_cat", "dv_pal_sleep_abn_cat", "dv_vdb_sleep_abn_cat", "dv_winkler_sleep_abn_cat")
We now convert variables into factors as appropriate, then specify
the imputation method as predictive mean matching (pmm) for
all variables.
# Merging sparse categories to aid convergence
data_merged %<>% mutate(
# Binarising e216a (mother's regular job since child's birth) to job vs no job
e216a = case_when(
e216a %in% c("F+Pt Job", "Ft Job", "Pt Job") ~ "Job",
e216a == "Other" | is.na(e216a) ~ NA_character_,
TRUE ~ e216a
),
# Combining classes I and II, and sending "other" categories to NA, in parental social class measures
fclrg90 = if_else(
fclrg90 %in% c("I Professional", "II Managerial"),
"I/II Professional/Managerial",
fclrg90
),
a0014 = case_when(
a0014 %in% c("SC 1", "SC 2") ~ "SC 1/2",
a0014 %in% c("Other", "Unsupported") ~ NA_character_,
TRUE ~ a0014
),
dv_sc_age_5 = case_when(
dv_sc_age_5 %in% c("I", "II") ~ "I/II",
dv_sc_age_5 == "Student-Volun Wk" ~ NA_character_,
TRUE ~ dv_sc_age_5
),
dv_sc_age_16 = case_when(
dv_sc_age_16 %in% c("I", "II") ~ "I/II",
dv_sc_age_16 %in% c("Dead", "Student") ~ NA_character_,
TRUE ~ dv_sc_age_16
),
# Binarising school type to comprehensive vs non-comprehensive
BDSTYPE = case_when(
BDSTYPE %in% c("Comprehensive", "Secondary Modern/Technical", "Scottish LEA") ~ "Comprehensive",
BDSTYPE %in% c("Grammar", "Independent Private", "Independent Special", "LEA Special") ~ "Non-comprehensive",
BDSTYPE == "Other" | is.na(BDSTYPE) ~ NA_character_
),
# Binarising ameni (lack of amenities) to any lack of amenities vs none
ameni = if_else(
ameni %in% c("1 occasion", "2 occasions"),
"At least 1 occasion",
ameni
),
# Binarising crowd (household overcrowding) to up to 1 person/room vs more than 1
crowd = if_else(
crowd %in% c("Over 1 to 1.5", "Over 1.5 to 2", "Over 2"),
"Over 1",
crowd
)
)
# Converting binary independent/auxiliary variables
for(var in c("a0230", "a0255", "d016a", "e216a", "BDSTYPE", "emosup", "vote97", "B9SCQ17", "B9SCQ19", "B9TEN", "ameni", "crowd", "divorce", "prmnh", "sepmumbcs", "dv_father_schl", "dv_med_5", "dv_med_10", "dv_incsource", "dv_alcohol", "dv_training", "dv_participation", "dv_organisations")){
data_merged[, var] %<>% as.factor()
}
# Converting unordered categorical independent/auxiliary variables
for(var in c("a0014", "dv_sc_age_5", "dv_sc_age_16", "dv_marital")){
data_merged[, var] %<>% as.factor()
}
# Converting ordered categorical independent/auxiliary variables
data_merged$a0014 %<>% factor(
ordered = TRUE,
levels = c("SC 1", "SC 2", "SC 3 NM", "SC 3 M", "SC 4", "SC 5")
)
data_merged$a0043b %<>% factor(
ordered = TRUE,
levels = c("Non Smoker", "Stopped Pre-Preg", "Stopped Dur-Preg", "Ctl Smokers 1 - 4", "Ctl Smokers 5 - 14", "Ctl Smokers >= 15")
)
data_merged$brfed %<>% factor(
ordered = TRUE,
levels = c("Breastfed for over 1 month", "Breastfed for under 1 month", "Not breastfed")
)
data_merged$resmove %<>% factor(
ordered = TRUE,
levels = c("None", "1-3", "4 or more")
)
data_merged$tenure %<>% factor(
ordered = TRUE,
levels = c("Owned at both time points", "Owned at one time point", "Rented at both time points")
)
data_merged$fclrg90 %<>% factor(
ordered = TRUE,
levels = c("I/II Professional/Managerial", "IIINM Skilled non-manual", "IIIM Skilled manual", "IV Partly skilled", "V Unskilled")
)
data_merged$dv_sc_age_5 %<>% factor(
ordered = TRUE,
levels = c("I/II", "IIINM", "IIIM", "IV", "V")
)
data_merged$dv_sc_age_16 %<>% factor(
ordered = TRUE,
levels = c("I/II", "III non-manual", "III manual", "IV", "V")
)
# Converting to factors, and making "normal" the reference category for, the multinomial regression variables
for(var in multinom_list){
data_merged[, var] %<>% as.factor() %>% relevel("Normal")
}
# Making a copy of the dataset, deleting the multinomial regression variables from the original version, and deleting the Poisson regression variables from the new version, so that they can be imputed separately
## The binary and ternary versions of each variable are highly collinear and therefore imputing them together causes chaos
data_merged_multi <- data_merged
data_merged %<>% select(-contains("abn_cat"))
data_merged_multi %<>% select(-(contains("sleep_abn") & !contains("abn_cat")))
# Creating vector of multiple imputation methods and setting method to pmm for all variables
method <- make.method(data_merged)
method[] <- "pmm"
method_multi <- make.method(data_merged_multi)
method_multi[] <- "pmm"
# Setting BCSID to not be imputed and removing it from the predictor set
method["bcsid"] <- ""
method_multi["bcsid"] <- ""
pred <- make.predictorMatrix(data_merged)
pred_multi <- make.predictorMatrix(data_merged_multi)
pred[, "bcsid"] <- 0
pred_multi[, "bcsid"] <- 0
We can now use the mice package to create multiple
imputed datasets from the data_merged data frame, which
contains missing data.
# Establishing a cluster for parallel computation to save time
## NB: THIS SECTION IS CURRENTLY SET TO CREATE AS MANY NODES AS THERE ARE CPU CORES AVAILABLE. THIS MAY CAUSE PROBLEMS IF THERE ARE OTHER TASKS RUNNING SIMULTANEOUSLY
### if you need to reduce the number, substitute `detectCores()` with e.g. `detectCores() - 2` or `6` in BOTH appearances below
cl <- makeCluster(detectCores())
# Setting random seed
clusterSetRNGStream(cl, 1858)
# Making the data frame visible in the cluster environment
clusterExport(cl, c("data_merged", "data_merged_multi", "method", "method_multi", "pred", "pred_multi"))
# Loading the mice package in the cluster environment
clusterEvalQ(cl, library(mice))
## [[1]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[2]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[3]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[4]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[5]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[6]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[7]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[8]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[9]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[10]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[11]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[12]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[13]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[14]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[15]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[16]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[17]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[18]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[19]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[20]]
## [1] "mice" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
# Running multiple imputation for Poisson variables dataset with 10 imputations per core used and a maximum of 15 iterations
## This code was run with 20 cores and therefore there are 200 imputations in total
imp_pars <-
parLapply(cl = cl, X = 1:detectCores(), fun = function(no){
mice(data_merged, vis = "monotone", method = method, predictorMatrix = pred, m = 10, maxit = 15, printFlag = FALSE)
})
# Merging imputations into a single mids object
imp_merged <- imp_pars[[1]]
for(i in 2:length(imp_pars)){
imp_merged <- mice::ibind(imp_merged, imp_pars[[i]])
}
rm(imp_pars)
# Running multiple imputation for multinomial variables dataset with 10 imputations per core used and a maximum of 15 iterations
## This code was run with 20 cores and therefore there are 200 imputations in total
imp_pars_multi <-
parLapply(cl = cl, X = 1:detectCores(), fun = function(no){
mice(data_merged_multi, vis = "monotone", method = method_multi, predictorMatrix = pred_multi, m = 10, maxit = 15, printFlag = FALSE)
})
stopCluster(cl)
# Merging imputations into a single mids object
imp_merged_multi <- imp_pars_multi[[1]]
for(i in 2:length(imp_pars_multi)){
imp_merged_multi <- mice::ibind(imp_merged_multi, imp_pars_multi[[i]])
}
rm(imp_pars_multi)
# Plotting trace lines to verify convergence
plot(imp_merged)
plot(imp_merged_multi)
We can now perform the regression analyses.
#######
# Function : pois_regr
# Input : ind_vars, a string containing a formula which describes the independent variables in the form "ind1 + ind2 + ind3 + ..." and var_of_interest, a string containing the variable of interest (e.g. "ind1")
# Method : The function fits modified Poisson regression models to both complete cases (for the age 5 and 10 exposures only) and the multiply imputed data using each of the dependent variables in turn and the specified independent variables, obtains the risk ratios, confidence intervals and p-values for the variable of interest from the model and binds these into a data frame
# Output : A data frame containing the risk ratios, confidence intervals and p-values obtained from modified Poisson regressions of the specified independent variable against the outcome variables, controlling for the other independent variables; these are derived from both complete cases and imputed data where possible, and the minimum E-value needed to explain away the association is also included
#######
pois_regr <- function(ind_vars, var_of_interest){
suppressWarnings({
poisson_regression <- function(dep_var, ind_vars, var_of_interest){
imp_model_fit <- with(imp_merged, coeftest(glm(as.formula(paste(dep_var, "~", ind_vars)), family = poisson(link = "log")), vcov. = sandwich))
imp_results <- summary(pool(imp_model_fit))
ind_var <- var_of_interest
RR <- exp(imp_results[imp_results$term == var_of_interest,]$estimate)
RR.LCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate + qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
RR.UCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate - qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
p <- imp_results[imp_results$term == var_of_interest,]$p.value
Evalue <- if(RR.LCI < 1 && RR.UCI > 1) 1 else if(RR.UCI < 1) (1/RR.UCI) + sqrt((1/RR.UCI)*((1/RR.UCI) - 1)) else if(RR.LCI > 1) RR.LCI + sqrt(RR.LCI*(RR.LCI - 1))
results_frame <- data.frame(dep_var, ind_var, RR, RR.LCI, RR.UCI, p, Evalue)
return(results_frame)
}
results <- lapply(pois_list, function(X) poisson_regression(X, ind_vars, var_of_interest)) %>% reduce(rbind)
return(results)
})
}
#######
# Function : multinom_regr
# Input : ind_vars, a string containing a formula which describes the independent variables in the form "ind1 + ind2 + ind3 + ..." (the dependent variable must be a factor) and var_of_interest, a string containing the variable of interest (e.g. "ind1")
# Method : The function fits multinomial logistic regression models to both complete cases (for the age 5 and 10 exposures only) and the multiply imputed data using each of the dependent variables in turn and the specified independent variables, obtains the relative risk ratios, confidence intervals and p-values for the variable of interest from the model and binds these into a data frame
# Output : A data frame containing the relative risk ratios, confidence intervals and p-values obtained from multinomial regressions of the specified independent variable against the dependent variable, controlling for the other independent variables; these are derived from both complete cases and imputed data where possible, and the minimum E-value needed to explain away the association is also included
#######
multinom_regr <- function(ind_vars, var_of_interest){
suppressWarnings({
multinomial_regression <- function(dep_var, ind_vars, var_of_interest){
imp_model_fit <- with(imp_merged_multi, multinom(as.formula(paste(dep_var, "~", ind_vars)), maxit = 200, trace = FALSE))
imp_results <- summary(pool(imp_model_fit))
ind_var <- var_of_interest
group <- imp_results[imp_results$term == var_of_interest,]$y.level
RRR <- exp(imp_results[imp_results$term == var_of_interest,]$estimate)
RRR.LCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate + qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
RRR.UCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate - qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
p <- imp_results[imp_results$term == var_of_interest,]$p.value
Evalue <- sapply(1:2, function(i){if(RRR.LCI[i] < 1 && RRR.UCI[i] > 1) 1 else if(RRR.UCI[i] < 1) (1/RRR.UCI[i]) + sqrt((1/RRR.UCI[i])*((1/RRR.UCI[i]) - 1)) else if(RRR.LCI[i] > 1) RRR.LCI[i] + sqrt(RRR.LCI[i]*(RRR.LCI[i] - 1))})
results_frame <- data.frame(dep_var, ind_var, group, RRR, RRR.LCI, RRR.UCI, p, Evalue)
return(results_frame)
}
results <- lapply(multinom_list, function(X) multinomial_regression(X, ind_vars, var_of_interest)) %>% reduce(rbind)
return(results)
})
}
# Setting random seed
set.seed(7917)
# Performing regressions against abnormal sleep
# Rutter (age 5)
rutter_5 <- pois_regr("d119 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + tenure + prmnh + crowd + ameni + brfed", "d119")
## Original model not retained as part of coeftest object. For additional model summary information (r.squared, df, etc.), consider passing `glance.coeftest()` an object where the underlying model has been saved, i.e.`lmtest::coeftest(..., save = TRUE)`.
## This message is displayed once per session.
print(rutter_5 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn d119 1.696 1.171 2.457 0.005 1.617
## 2 dv_sd_sleep_abn d119 1.533 0.939 2.503 0.089 1.000
## 3 dv_pal_sleep_abn d119 0.881 0.677 1.145 0.343 1.000
## 4 dv_vdb_sleep_abn d119 1.203 0.879 1.645 0.249 1.000
## 5 dv_winkler_sleep_abn d119 1.784 1.296 2.455 0.000 1.915
rutter_5_multi <- multinom_regr("d119 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + tenure + prmnh + crowd + ameni + brfed", "d119")
print(rutter_5_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var group RRR RRR.LCI RRR.UCI p Evalue
## 1 dv_sr_sleep_abn_cat d119 Long 1.798 0.526 6.146 0.349 1.000
## 2 dv_sr_sleep_abn_cat d119 Short 1.891 1.180 3.030 0.008 1.642
## 3 dv_sd_sleep_abn_cat d119 Long 1.543 0.708 3.363 0.276 1.000
## 4 dv_sd_sleep_abn_cat d119 Short 1.547 0.847 2.825 0.157 1.000
## 5 dv_pal_sleep_abn_cat d119 Long 0.854 0.568 1.282 0.446 1.000
## 6 dv_pal_sleep_abn_cat d119 Short 0.878 0.473 1.632 0.681 1.000
## 7 dv_vdb_sleep_abn_cat d119 Long 1.250 0.770 2.029 0.367 1.000
## 8 dv_vdb_sleep_abn_cat d119 Short 1.363 0.565 3.289 0.491 1.000
## 9 dv_winkler_sleep_abn_cat d119 Long 1.172 0.572 2.400 0.665 1.000
## 10 dv_winkler_sleep_abn_cat d119 Short 2.352 1.385 3.995 0.002 2.115
# Rutter (age 10)
rutter_10 <- pois_regr("BD3MRUTT + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD3MRUTT")
print(rutter_10 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn BD3MRUTT 1.688 1.048 2.718 0.032 1.271
## 2 dv_sd_sleep_abn BD3MRUTT 2.115 1.194 3.748 0.011 1.675
## 3 dv_pal_sleep_abn BD3MRUTT 1.016 0.741 1.392 0.923 1.000
## 4 dv_vdb_sleep_abn BD3MRUTT 1.324 0.898 1.952 0.158 1.000
## 5 dv_winkler_sleep_abn BD3MRUTT 1.313 0.913 1.888 0.143 1.000
rutter_10_multi <- multinom_regr("BD3MRUTT + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD3MRUTT")
print(rutter_10_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var group RRR RRR.LCI RRR.UCI p Evalue
## 1 dv_sr_sleep_abn_cat BD3MRUTT Long 0.927 0.190 4.536 0.926 1.000
## 2 dv_sr_sleep_abn_cat BD3MRUTT Short 2.139 1.235 3.706 0.007 1.773
## 3 dv_sd_sleep_abn_cat BD3MRUTT Long 1.853 0.709 4.841 0.208 1.000
## 4 dv_sd_sleep_abn_cat BD3MRUTT Short 2.446 1.181 5.065 0.017 1.643
## 5 dv_pal_sleep_abn_cat BD3MRUTT Long 1.023 0.627 1.669 0.928 1.000
## 6 dv_pal_sleep_abn_cat BD3MRUTT Short 1.149 0.538 2.457 0.720 1.000
## 7 dv_vdb_sleep_abn_cat BD3MRUTT Long 1.619 0.852 3.076 0.143 1.000
## 8 dv_vdb_sleep_abn_cat BD3MRUTT Short 1.904 0.652 5.559 0.239 1.000
## 9 dv_winkler_sleep_abn_cat BD3MRUTT Long 1.034 0.447 2.395 0.938 1.000
## 10 dv_winkler_sleep_abn_cat BD3MRUTT Short 1.420 0.765 2.635 0.268 1.000
# Child Development Scale
cds <- pois_regr("dv_cds_10 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "dv_cds_10")
print(cds %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn dv_cds_10 3.578 2.042 6.267 0.000 3.501
## 2 dv_sd_sleep_abn dv_cds_10 3.972 1.991 7.925 0.000 3.395
## 3 dv_pal_sleep_abn dv_cds_10 0.822 0.562 1.204 0.315 1.000
## 4 dv_vdb_sleep_abn dv_cds_10 1.114 0.741 1.676 0.604 1.000
## 5 dv_winkler_sleep_abn dv_cds_10 2.304 1.461 3.635 0.000 2.281
cds_multi <- multinom_regr("dv_cds_10 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "dv_cds_10")
print(cds_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var group RRR RRR.LCI RRR.UCI p Evalue
## 1 dv_sr_sleep_abn_cat dv_cds_10 Long 2.224 0.438 11.296 0.335 1.000
## 2 dv_sr_sleep_abn_cat dv_cds_10 Short 5.071 2.628 9.784 0.000 4.696
## 3 dv_sd_sleep_abn_cat dv_cds_10 Long 2.734 0.995 7.517 0.052 1.000
## 4 dv_sd_sleep_abn_cat dv_cds_10 Short 4.568 2.077 10.049 0.000 3.572
## 5 dv_pal_sleep_abn_cat dv_cds_10 Long 0.855 0.465 1.572 0.614 1.000
## 6 dv_pal_sleep_abn_cat dv_cds_10 Short 1.055 0.401 2.775 0.914 1.000
## 7 dv_vdb_sleep_abn_cat dv_cds_10 Long 1.279 0.619 2.643 0.506 1.000
## 8 dv_vdb_sleep_abn_cat dv_cds_10 Short 1.610 0.515 5.031 0.413 1.000
## 9 dv_winkler_sleep_abn_cat dv_cds_10 Long 0.985 0.381 2.547 0.975 1.000
## 10 dv_winkler_sleep_abn_cat dv_cds_10 Short 3.727 1.730 8.029 0.001 2.854
# Malaise Inventory
malaise <- pois_regr("BD4MAL + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD4MAL")
print(malaise %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn BD4MAL 3.602 1.839 7.055 0.000 3.081
## 2 dv_sd_sleep_abn BD4MAL 1.575 0.818 3.034 0.176 1.000
## 3 dv_pal_sleep_abn BD4MAL 1.158 0.787 1.703 0.457 1.000
## 4 dv_vdb_sleep_abn BD4MAL 1.178 0.761 1.825 0.462 1.000
## 5 dv_winkler_sleep_abn BD4MAL 1.492 0.964 2.311 0.074 1.000
malaise_multi <- multinom_regr("BD4MAL + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD4MAL")
print(malaise_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var group RRR RRR.LCI RRR.UCI p Evalue
## 1 dv_sr_sleep_abn_cat BD4MAL Long 2.219 0.203 24.200 0.514 1.000
## 2 dv_sr_sleep_abn_cat BD4MAL Short 5.554 2.461 12.536 0.000 4.357
## 3 dv_sd_sleep_abn_cat BD4MAL Long 1.492 0.321 6.937 0.610 1.000
## 4 dv_sd_sleep_abn_cat BD4MAL Short 1.622 0.611 4.307 0.332 1.000
## 5 dv_pal_sleep_abn_cat BD4MAL Long 1.255 0.722 2.180 0.422 1.000
## 6 dv_pal_sleep_abn_cat BD4MAL Short 1.367 0.551 3.393 0.501 1.000
## 7 dv_vdb_sleep_abn_cat BD4MAL Long 1.165 0.577 2.354 0.670 1.000
## 8 dv_vdb_sleep_abn_cat BD4MAL Short 1.114 0.325 3.820 0.864 1.000
## 9 dv_winkler_sleep_abn_cat BD4MAL Long 1.482 0.663 3.311 0.338 1.000
## 10 dv_winkler_sleep_abn_cat BD4MAL Short 1.624 0.802 3.288 0.179 1.000
# Behavioural or emotional problems
beh_emo <- pois_regr("rd6m_1 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "rd6m_1")
print(beh_emo %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn rd6m_1 1.612 1.113 2.335 0.012 1.467
## 2 dv_sd_sleep_abn rd6m_1 1.518 0.933 2.471 0.094 1.000
## 3 dv_pal_sleep_abn rd6m_1 1.148 0.846 1.560 0.376 1.000
## 4 dv_vdb_sleep_abn rd6m_1 0.909 0.621 1.331 0.626 1.000
## 5 dv_winkler_sleep_abn rd6m_1 1.049 0.739 1.488 0.790 1.000
beh_emo_multi <- multinom_regr("rd6m_1 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "rd6m_1")
print(beh_emo_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var group RRR RRR.LCI RRR.UCI p Evalue
## 1 dv_sr_sleep_abn_cat rd6m_1 Long 1.312 0.259 6.655 0.743 1.000
## 2 dv_sr_sleep_abn_cat rd6m_1 Short 1.852 1.076 3.188 0.027 1.363
## 3 dv_sd_sleep_abn_cat rd6m_1 Long 1.149 0.489 2.699 0.750 1.000
## 4 dv_sd_sleep_abn_cat rd6m_1 Short 2.005 0.984 4.085 0.057 1.000
## 5 dv_pal_sleep_abn_cat rd6m_1 Long 1.451 0.805 2.614 0.217 1.000
## 6 dv_pal_sleep_abn_cat rd6m_1 Short 1.917 0.828 4.439 0.130 1.000
## 7 dv_vdb_sleep_abn_cat rd6m_1 Long 1.016 0.494 2.091 0.965 1.000
## 8 dv_vdb_sleep_abn_cat rd6m_1 Short 1.074 0.402 2.869 0.887 1.000
## 9 dv_winkler_sleep_abn_cat rd6m_1 Long 0.916 0.412 2.037 0.830 1.000
## 10 dv_winkler_sleep_abn_cat rd6m_1 Short 1.080 0.576 2.026 0.810 1.000
# Internalising behaviour in childhood
int_beh <- pois_regr("intbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "intbcsz")
print(int_beh %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn intbcsz 1.696 1.051 2.737 0.031 1.284
## 2 dv_sd_sleep_abn intbcsz 1.508 0.912 2.493 0.110 1.000
## 3 dv_pal_sleep_abn intbcsz 1.129 0.822 1.550 0.456 1.000
## 4 dv_vdb_sleep_abn intbcsz 1.212 0.864 1.699 0.266 1.000
## 5 dv_winkler_sleep_abn intbcsz 1.189 0.848 1.668 0.316 1.000
int_beh_multi <- multinom_regr("intbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "intbcsz")
print(int_beh_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var group RRR RRR.LCI RRR.UCI p Evalue
## 1 dv_sr_sleep_abn_cat intbcsz Long 2.350 0.416 13.278 0.334 1.000
## 2 dv_sr_sleep_abn_cat intbcsz Short 1.930 1.080 3.448 0.027 1.375
## 3 dv_sd_sleep_abn_cat intbcsz Long 1.134 0.489 2.633 0.769 1.000
## 4 dv_sd_sleep_abn_cat intbcsz Short 1.903 0.965 3.756 0.064 1.000
## 5 dv_pal_sleep_abn_cat intbcsz Long 1.259 0.791 2.005 0.332 1.000
## 6 dv_pal_sleep_abn_cat intbcsz Short 1.481 0.690 3.177 0.314 1.000
## 7 dv_vdb_sleep_abn_cat intbcsz Long 1.429 0.768 2.659 0.261 1.000
## 8 dv_vdb_sleep_abn_cat intbcsz Short 1.473 0.557 3.892 0.435 1.000
## 9 dv_winkler_sleep_abn_cat intbcsz Long 1.096 0.516 2.330 0.811 1.000
## 10 dv_winkler_sleep_abn_cat intbcsz Short 1.298 0.706 2.386 0.401 1.000
# Externalising behaviour in childhood
ext_beh <- pois_regr("extbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "extbcsz")
print(ext_beh %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn extbcsz 2.280 1.461 3.557 0.000 2.282
## 2 dv_sd_sleep_abn extbcsz 1.669 0.948 2.937 0.077 1.000
## 3 dv_pal_sleep_abn extbcsz 1.053 0.734 1.511 0.778 1.000
## 4 dv_vdb_sleep_abn extbcsz 1.336 0.879 2.029 0.176 1.000
## 5 dv_winkler_sleep_abn extbcsz 1.779 1.248 2.536 0.002 1.805
ext_beh_multi <- multinom_regr("extbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "extbcsz")
print(ext_beh_multi %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var group RRR RRR.LCI RRR.UCI p Evalue
## 1 dv_sr_sleep_abn_cat extbcsz Long 4.837 1.151 20.334 0.032 1.567
## 2 dv_sr_sleep_abn_cat extbcsz Short 2.821 1.558 5.107 0.001 2.491
## 3 dv_sd_sleep_abn_cat extbcsz Long 1.194 0.495 2.880 0.693 1.000
## 4 dv_sd_sleep_abn_cat extbcsz Short 2.013 0.958 4.232 0.066 1.000
## 5 dv_pal_sleep_abn_cat extbcsz Long 1.128 0.620 2.053 0.693 1.000
## 6 dv_pal_sleep_abn_cat extbcsz Short 1.352 0.534 3.422 0.526 1.000
## 7 dv_vdb_sleep_abn_cat extbcsz Long 1.666 0.784 3.543 0.186 1.000
## 8 dv_vdb_sleep_abn_cat extbcsz Short 1.914 0.690 5.309 0.213 1.000
## 9 dv_winkler_sleep_abn_cat extbcsz Long 1.126 0.513 2.474 0.767 1.000
## 10 dv_winkler_sleep_abn_cat extbcsz Short 2.452 1.173 5.128 0.018 1.623
We now write the results of the regression analyses to .csv files.
fwrite(rutter_5, ".\\Results\\rutter_5.csv")
fwrite(rutter_10, ".\\Results\\rutter_10.csv")
fwrite(cds, ".\\Results\\cds.csv")
fwrite(malaise, ".\\Results\\malaise.csv")
fwrite(beh_emo, ".\\Results\\beh_emo.csv")
fwrite(int_beh, ".\\Results\\int_beh.csv")
fwrite(ext_beh, ".\\Results\\ext_beh.csv")
fwrite(rutter_5_multi, ".\\Results\\rutter_5_multi.csv")
fwrite(rutter_10_multi, ".\\Results\\rutter_10_multi.csv")
fwrite(cds_multi, ".\\Results\\cds_multi.csv")
fwrite(malaise_multi, ".\\Results\\malaise_multi.csv")
fwrite(beh_emo_multi, ".\\Results\\beh_emo_multi.csv")
fwrite(int_beh_multi, ".\\Results\\int_beh_multi.csv")
fwrite(ext_beh_multi, ".\\Results\\ext_beh_multi.csv")
We suspect that the association between early-life mental health and sleep in adulthood may be largely mediated by adult mental health. To test this, we informally perform mediation path analyses. We have already established that all seven of the exposures are associated with at least some of the adult sleep variables, but we still need to determine (a) whether childhood mental health is associated with adult mental health (here we are using Malaise and WEMWBS scores at age 42) and (b) whether controlling for these adult mental health variables attenuates the associations we found among childhood mental health and adult sleep variables.
#######
# Function : med_pois_regr
# Input : ind_vars, a string containing a formula which describes the independent variables in the form "ind1 + ind2 + ind3 + ..." and var_of_interest, a string containing the variable of interest (e.g. "ind1")
# Method : The function fits modified Poisson regression models to the multiply imputed data using each of the mediating adult mental health variables in turn and the independent variables, obtains the risk ratios, confidence intervals and p-values for the variable of interest from the model and binds these into a data frame
# Output : A data frame containing the risk ratios, confidence intervals and p-values obtained from modified Poisson regressions of the mediating adult mental health variables against the independent variables, controlling for the other independent variables
#######
med_pois_regr <- function(ind_vars, var_of_interest){
suppressWarnings({
mediator_poisson_regression <- function(dep_var, ind_vars, var_of_interest){
imp_model_fit <- with(imp_merged, coeftest(glm(as.formula(paste(dep_var, "~", ind_vars)), family = poisson(link = "log")), vcov. = sandwich))
imp_results <- summary(pool(imp_model_fit))
ind_var <- var_of_interest
impRR <- exp(imp_results[imp_results$term == var_of_interest,]$estimate)
impRR.LCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate + qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
impRR.UCI <- exp(imp_results[imp_results$term == var_of_interest,]$estimate - qnorm(0.05/2)*imp_results[imp_results$term == var_of_interest,]$std.error)
impP <- imp_results[imp_results$term == var_of_interest,]$p.value
results_frame <- data.frame(dep_var, ind_var, impRR, impRR.LCI, impRR.UCI, impP)
return(results_frame)
}
results <- lapply(c("BD9MAL", "BD9WEMWB"), function(X) mediator_poisson_regression(X, ind_vars, var_of_interest)) %>% reduce(rbind)
return(results)
})
}
# Rutter (age 5)
rutter_5_med1 <- med_pois_regr("d119 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + tenure + prmnh + crowd + ameni + brfed", "d119")
rutter_5_med2 <- pois_regr("d119 + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + tenure + prmnh + crowd + ameni + brfed", "d119")
print(rutter_5_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1 BD9MAL d119 1.818 1.537 2.151 0
## 2 BD9WEMWB d119 0.837 0.802 0.875 0
print(rutter_5_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn d119 1.264 0.881 1.813 0.204 1.000
## 2 dv_sd_sleep_abn d119 1.276 0.788 2.065 0.323 1.000
## 3 dv_pal_sleep_abn d119 0.869 0.667 1.133 0.301 1.000
## 4 dv_vdb_sleep_abn d119 1.149 0.841 1.571 0.384 1.000
## 5 dv_winkler_sleep_abn d119 1.680 1.226 2.303 0.001 1.753
# Rutter (age 10)
rutter_10_med1 <- med_pois_regr("BD3MRUTT + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD3MRUTT")
rutter_10_med2 <- pois_regr("BD3MRUTT + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD3MRUTT")
print(rutter_10_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1 BD9MAL BD3MRUTT 2.317 1.886 2.847 0
## 2 BD9WEMWB BD3MRUTT 0.764 0.724 0.807 0
print(rutter_10_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn BD3MRUTT 1.103 0.696 1.748 0.676 1
## 2 dv_sd_sleep_abn BD3MRUTT 1.640 0.940 2.861 0.082 1
## 3 dv_pal_sleep_abn BD3MRUTT 1.001 0.731 1.372 0.994 1
## 4 dv_vdb_sleep_abn BD3MRUTT 1.247 0.841 1.848 0.273 1
## 5 dv_winkler_sleep_abn BD3MRUTT 1.208 0.839 1.737 0.310 1
# Child Development Scale
cds_med1 <- med_pois_regr("dv_cds_10 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "dv_cds_10")
cds_med2 <- pois_regr("dv_cds_10 + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + bmi + fclrg90 + tenure + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "dv_cds_10")
print(cds_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1 BD9MAL dv_cds_10 1.777 1.390 2.273 0
## 2 BD9WEMWB dv_cds_10 0.847 0.794 0.903 0
print(cds_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn dv_cds_10 2.680 1.551 4.631 0.000 2.475
## 2 dv_sd_sleep_abn dv_cds_10 3.355 1.718 6.552 0.000 2.829
## 3 dv_pal_sleep_abn dv_cds_10 0.813 0.554 1.191 0.289 1.000
## 4 dv_vdb_sleep_abn dv_cds_10 1.068 0.713 1.600 0.750 1.000
## 5 dv_winkler_sleep_abn dv_cds_10 2.189 1.389 3.449 0.001 2.124
# Malaise Inventory
malaise_med1 <- med_pois_regr("BD4MAL + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD4MAL")
malaise_med2 <- pois_regr("BD4MAL + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "BD4MAL")
print(malaise_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1 BD9MAL BD4MAL 7.002 3.986 12.302 0
## 2 BD9WEMWB BD4MAL 0.668 0.588 0.759 0
print(malaise_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn BD4MAL 1.677 0.890 3.162 0.111 1
## 2 dv_sd_sleep_abn BD4MAL 0.997 0.505 1.967 0.993 1
## 3 dv_pal_sleep_abn BD4MAL 1.130 0.747 1.708 0.563 1
## 4 dv_vdb_sleep_abn BD4MAL 1.030 0.652 1.628 0.898 1
## 5 dv_winkler_sleep_abn BD4MAL 1.270 0.799 2.019 0.314 1
# Behavioural or emotional problems
beh_emo_med1 <- med_pois_regr("rd6m_1 + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "rd6m_1")
beh_emo_med2 <- pois_regr("rd6m_1 + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "rd6m_1")
print(beh_emo_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1 BD9MAL rd6m_1 1.404 1.159 1.701 0.001
## 2 BD9WEMWB rd6m_1 0.906 0.848 0.968 0.004
print(beh_emo_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn rd6m_1 1.350 0.950 1.918 0.095 1
## 2 dv_sd_sleep_abn rd6m_1 1.372 0.846 2.226 0.201 1
## 3 dv_pal_sleep_abn rd6m_1 1.142 0.840 1.551 0.398 1
## 4 dv_vdb_sleep_abn rd6m_1 0.884 0.604 1.294 0.527 1
## 5 dv_winkler_sleep_abn rd6m_1 1.011 0.713 1.433 0.952 1
# Internalising behaviour in childhood
int_beh_med1 <- med_pois_regr("intbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "intbcsz")
int_beh_med2 <- pois_regr("intbcsz + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "intbcsz")
print(int_beh_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1 BD9MAL intbcsz 3.451 2.863 4.160 0
## 2 BD9WEMWB intbcsz 0.743 0.705 0.784 0
print(int_beh_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn intbcsz 0.978 0.612 1.564 0.927 1
## 2 dv_sd_sleep_abn intbcsz 1.111 0.667 1.848 0.687 1
## 3 dv_pal_sleep_abn intbcsz 1.109 0.804 1.530 0.527 1
## 4 dv_vdb_sleep_abn intbcsz 1.116 0.785 1.586 0.540 1
## 5 dv_winkler_sleep_abn intbcsz 1.061 0.753 1.493 0.736 1
# Externalising behaviour in childhood
ext_beh_med1 <- med_pois_regr("extbcsz + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "extbcsz")
ext_beh_med2 <- pois_regr("extbcsz + BD9MAL + BD9WEMWB + a0005a + a0043b + a0195a + a0278 + a0014 + dv_par_edu_birth + dv_sc_age_5 + e216a + dv_cog_abil_5 + dv_med_5 + d016a + dv_bas_g + dv_med_10 + dv_sc_age_16 + bmi + fclrg90 + tenure + divorce + sepmumbcs + prmnh + crowd + ameni + brfed + resmove", "extbcsz")
print(ext_beh_med1 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var impRR impRR.LCI impRR.UCI impP
## 1 BD9MAL extbcsz 2.435 2.023 2.932 0
## 2 BD9WEMWB extbcsz 0.759 0.709 0.812 0
print(ext_beh_med2 %>% mutate_if(is.numeric, function(x) round(x, 3)))
## dep_var ind_var RR RR.LCI RR.UCI p Evalue
## 1 dv_sr_sleep_abn extbcsz 1.505 0.973 2.327 0.067 1.000
## 2 dv_sd_sleep_abn extbcsz 1.300 0.731 2.312 0.373 1.000
## 3 dv_pal_sleep_abn extbcsz 1.038 0.720 1.495 0.843 1.000
## 4 dv_vdb_sleep_abn extbcsz 1.255 0.824 1.909 0.291 1.000
## 5 dv_winkler_sleep_abn extbcsz 1.644 1.155 2.340 0.006 1.578
# Exporting mediation analysis results
fwrite(rutter_5_med1, ".\\Results\\Mediation\\rutter_5_med1.csv")
fwrite(rutter_5_med2, ".\\Results\\Mediation\\rutter_5_med2.csv")
fwrite(rutter_10_med1, ".\\Results\\Mediation\\rutter_10_med1.csv")
fwrite(rutter_10_med2, ".\\Results\\Mediation\\rutter_10_med2.csv")
fwrite(cds_med1, ".\\Results\\Mediation\\cds_med1.csv")
fwrite(cds_med2, ".\\Results\\Mediation\\cds_med2.csv")
fwrite(malaise_med1, ".\\Results\\Mediation\\malaise_med1.csv")
fwrite(malaise_med2, ".\\Results\\Mediation\\malaise_med2.csv")
fwrite(beh_emo_med1, ".\\Results\\Mediation\\beh_emo_med1.csv")
fwrite(beh_emo_med2, ".\\Results\\Mediation\\beh_emo_med2.csv")
fwrite(int_beh_med1, ".\\Results\\Mediation\\int_beh_med1.csv")
fwrite(int_beh_med2, ".\\Results\\Mediation\\int_beh_med2.csv")
fwrite(ext_beh_med1, ".\\Results\\Mediation\\ext_beh_med1.csv")
fwrite(ext_beh_med2, ".\\Results\\Mediation\\ext_beh_med2.csv")